In [ ]:
###############
#NOTES
###############
"""
I need to look into filtering functions, obviously some of the data is fucked up,
for example sat 5 viewed from mah3 is mostly normal data but then a few values are
10^9 so whats up with that
FILTERING
as for receiver bias determination, the paper actually says that for sites above 65 degrees
shouldn't use minimum scalloping and should use zero tec method and least squares instead
the problem is that it would be better if I had night time data, my data is only from 6am to
about 1pm
next week I guess I should look into least squares receiver bias estimation and zero tec
or i could combine least squares and minimum scalloping like they do in the paper
my shit needs to be faster
"""
In [2]:
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from matplotlib.dates import YearLocator, MonthLocator, DateFormatter
from pandas import DataFrame, Panel4D, read_hdf
from glob import glob
from datetime import datetime
import time
import matplotlib
import random
import gps
from IPython.display import Image
from pandas import DataFrame,Series
import CoordTransforms3 as CT
In [3]:
data = read_hdf('svn23.h5','data')
start = datetime(year=2015,month=10,day=7,hour=6,minute=19,second=0)
end = datetime(year=2015,month=10,day=7,hour=6,minute=23,second=0)
plt.figure(figsize=(16,16))
fmt = DateFormatter('%H:%M:%S')
plt.subplot(311).xaxis.set_major_formatter(fmt)
colors = ['b','r','g','m','c','y']
extra = [a for a in matplotlib.colors.cnames]
biasobj = gps.satelliteBias('jplg2800.15i','P1C11510.DCB',None)
bias = biasobj.dict[(23,1)]
sites = list(data.labels)
sites.remove('mah72800.15o')
sites.remove('mah92800.15o')
labs=[]
for site in sites:
try:
color=colors.pop()
except:
color=random.choice(extra)
ranges = gps.getRangesalt(data,site,maxjump=1.0)
plots=[]
for drange in ranges:
tec,err = gps.getTecalt(data,site,drange)
tec-=bias
plots.append(tec)
for p in plots:
if site[:4] not in labs:
plt.plot(p,'.',c=color,label=site[:4])
labs.append(site[:4])
else:
plt.plot(p,'.',c=color,label='')
plt.xlim([start,end])
plt.legend()
plt.grid()
plt.xlabel('time')
plt.ylabel('TECu')
plt.title('Trying to Replicate the Rideout Plot')
plt.show()
Image("gradient.jpg")
Out[3]:
This plot was made using most of the same methods from tec.py, I still need to convert to vTEC and maybe look into filtering? Why is there data missing from my plot? The data from mah6 is negative and then dissapears. Other lines are too high. This could be because I am plotting slant TEC, it could also be because I haven't accounted for receiver bias, It could also be because I didn't use the same constants. His data looks steeper though, look at mah3, our plots generally have the same shape but the spike in his data is more pronounced.
In [6]:
from mpl_toolkits.mplot3d import Axes3D
files = glob('/home/greg/Documents/Summer Research/rinex files/mah*')
head,data = gps.rinexobs(files[0],returnHead=True)
nav = gps.readRinexNav('./brdc2800.15n')
import datetime
base = datetime.datetime(2015,10,7,0)
arr = np.array([base + datetime.timedelta(minutes=i) for i in range(60*24)])
satpos1 = gps.getSatXYZ(nav,1,arr)
satpos2 = gps.getSatXYZ(nav,2,arr)
satpos3 = gps.getSatXYZ(nav,3,arr)
recpos=np.asarray(head['APPROX POSITION XYZ'])[:,None]
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111, projection='3d')
ax.set_aspect('equal')
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
dist = np.linalg.norm(recpos)
x = dist * np.outer(np.cos(u), np.sin(v))
y = dist * np.outer(np.sin(u), np.sin(v))
z = dist * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x,y,z,rstride=4,cstride=4,color='b',linewidth=1,alpha=1)
ax.plot(recpos[0],recpos[1],recpos[2],'r.',markersize=10)
ax.plot(satpos1[:,0],satpos1[:,1],satpos1[:,2],label='sv1')
ax.plot(satpos2[:,0],satpos2[:,1],satpos2[:,2],label='sv2')
ax.plot(satpos3[:,0],satpos3[:,1],satpos3[:,2],label='sv3')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()
In [73]:
start = datetime(2015,10,7,6,19,0)
end = datetime(2015,10,7,6,23,0)
data = read_hdf('svn23.h5','data')
nav = gps.readRinexNav('./brdc2800.15n')
recpos = np.array([[-2265800.1633],
[-1401757.874 ],
[ 5775924.6396]])
print('data loaded')
satpos = gps.getSatXYZ(nav,23,data.items)
print('sat positions in XYZ')
recposwgs = CT.ecef2wgs(recpos)
satvecenu = CT.ecef2enul(np.asarray(satpos,float),recposwgs.T)
satveccart = CT.enu2cartisian(satvecenu)
satvecsph = CT.cartisian2Sphereical(satveccart)
el = Series(satvecsph[2,:],index=data.items)
z = Series(gps.getZ(el),index=data.items)
print('mapping function')
plt.figure(figsize=(16,16))
fmt = DateFormatter('%H:%M:%S')
plt.subplot(311).xaxis.set_major_formatter(fmt)
colors = ['b','r','g','m','c','y']
biasobj = gps.satelliteBias('jplg2800.15i','P1C11510.DCB',None)
bias = biasobj.dict[(23,1)]
sites = list(data.labels)
sites.remove('mah72800.15o')
sites.remove('mah92800.15o')
for site in sites:
print(site[:4])
color=colors.pop()
ranges = gps.getRangesalt(data,site,maxjump=1.0)
tecs=[]
times=[]
for drange in ranges:
tec,err = gps.getTecalt(data,site,drange)
tec-=bias
for d in tec.index:
if(d in z.index):
tec[d]/=z[d]
tecs.append(tec)
times.append(tec.index)
vtec = np.hstack((p for p in tecs))
times=np.hstack((j for j in times))
vTEC = Series(vtec,index=times)
rbias = gps.minScalErr(vTEC,el)
print('rbias = ',rbias)
print('times = ',len(vTEC))
plt.plot(vTEC,'o-',c=color,label=site[:4],fillstyle='none')
plt.xlim([start,end])
plt.ylim([-5,30])
plt.legend()
plt.grid()
plt.xlabel('time')
plt.ylabel('VTEC in TECu')
plt.title('Replica of the Rideout Plot')
plt.show()
Image("gradient.jpg")
Out[73]:
In [7]:
#READ FILES
files = glob("/home/greg/Documents/Summer Research/rinex files/mah*")
head,data = gps.rinexobs(files[0],returnHead=True)
nav = gps.readRinexNav('./brdc2800.15n')
biasobj = gps.satelliteBias('jplg2800.15i','P1C11510.DCB',None)
In [8]:
bias = gps.minScalBias(head,data,nav,biasobj)
In [11]:
####################
#LEAST SQUARES
####################
sv=23
ranges = gps.getRanges(data,sv)
tecs=[]
times=[]
for drange in ranges:
tec,err = gps.getTec(data,sv,drange)
tec-=bias
tecs.append(tec)
times.append(tec.index)
stec = np.hstack((p for p in tecs))
times=np.hstack((j for j in times))
stec = Series(stec,index=times)
In [ ]: